function [log_likelihood, nan_pct, sim_factors, sim_l_contrib_all_factors, theta_i, kappa_i, p_type, sim_l_contrib_1j]=ml_risk_factors_score_types_cdf(theta)
%%% Joint optimization step for likelihood of observing measures and
%%% choices on lottery tasks



 global f_ijs X_ind_s measures ind_id factor_n  ind_n draw_n c_char_n b_char_n m_char_n cont_char_n thresh  X_n theta_xy lottery_n lottery_choice_ind n_type maxima minima sig_th_scale sig_th_lower_limit 
    
    % each measure is associated with a unique mean and loading coefficient
    % for the factor 
    f_measures=measures;
    columns=sum(c_char_n);
    % # of observations including rows with simulated factor draws for each
    % person
    rows=size(f_measures,1);
    
    sex=X_ind_s(:,1);
     
    %%% once lottery choices are added in, everything is jointly estimated,
    %%% so now for each simulated draw of all factors combined, need to
    %%% multiply the likelihoods together 
    sim_l_contrib_all_factors=NaN(rows, columns);
    
    %%% simulated likelihood contributions for next step (step2)
    sim_l_contrib_1j=NaN(rows,3);
    
    %%% stores the value of each factor based on Xalpha from characteristics
    %%% and corresponding factor draws
    sim_factors=zeros(rows,factor_n);
    
    %%% tracks the first column into which the simulated likelihoods for a
    %%% particular factor (1,2,...,n) are put in
    z=1;
    


    
    for i=1:factor_n
        
        % gets coeffs on the demographic vars for each factor 
        alphas=theta(1:X_n);
        theta(1:X_n)=[];


        % all multi-valued variables have the same amount of categories
        % for a given factor
        % number of nuisance parameters (thresh-mean) to be estimated for
        % this particular factor
        size_const=b_char_n(i)+thresh(i)*m_char_n(i)+cont_char_n(i);
        %%% intercept for 1st measure of each factor is normalized to 0
        constants=[0; theta(1:size_const-1)];
        %%% factor-specific intercept
        f_mean=theta(size_const);
        % one loading per measure
        loading=[1 theta(size_const+1:size_const+c_char_n(i)-1)'];
        %%% factor orthogonal component sd
        sd=exp(theta(size_const+c_char_n(i)));
        %%% error sd for the continous measures
        sd_cont=exp(theta(size_const+c_char_n(i)+1:size_const+c_char_n(i)+cont_char_n(i)));
        theta(1:size_const+c_char_n(i)+cont_char_n(i))=[];        
        %%% to be filled in gradually with the binary, multi, and
        %%% continuous measure likelihood contributions
        sim_l_contrib=NaN(rows,c_char_n(i));

        %%% "f_ij" are simulated factor draws (common to all measures)
        %%% associated with a given factor
        %%% gets the simulated value of factor_i 
        sim_factor=f_mean+X_ind_s*alphas+f_ijs(:,i)*sd;


        if b_char_n(i) > 0

            % gets the binary measures for the factor
            X_c_b=f_measures(:,1:b_char_n(i));
            %nuisance parameters for the binary measures associated with the
            %factor
            t_k_mean_b=constants(1:b_char_n(i))';
            % remove the already assigned parameters from the constant list
            constants(1:b_char_n(i))=[];
            % loadings for binary measure
            loading_b=loading(1:b_char_n(i));
            % remove the already assigned parameters from the loading list
            loading(1:b_char_n(i))=[];


            %%%%%%%%%%%% Simulated likelihood contributions of the BINARY
            %%%%%%%%%%%% measures (calculates sim_u and sim_l_contrib only for those now)

            %%% simulated measure value for each individual and simulated factor
            %%% draw - excludes the constant term for multi-valued categorical
            %%% estimation
            sim_u=sim_factor*loading_b;

            temp=bsxfun(@plus,-sim_u,t_k_mean_b);       
            sim_l_contrib(:,1:b_char_n(i))=((1-normcdf(temp)).^(X_c_b)).*(normcdf(temp).^(1-X_c_b));

        end
        
        
        
        
        
        
        %%%%%%%%%%%% Simulated likelihood contributions of the MULTI-VALUED
        %%%%%%%%%%%% measures (calculates sim_u and sim_l_contrib for those now)
        

        
        if m_char_n(i) > 0
            %%% gets the multi-valued measures for the factor
            X_c_m=f_measures(:,b_char_n(i)+1:b_char_n(i)+m_char_n(i));

            t_k_mean_raw=zeros(thresh(i),m_char_n(i));
            for j=1:thresh(i)
                t_k_mean_raw(j,:)=constants(1+(j-1)*m_char_n(i):j*m_char_n(i));
            end

            t_k_mean_m=zeros(thresh(i),m_char_n(i));
            t_k_mean_m(1,:)=t_k_mean_raw(1,:);
            % ensures monotonicity in threshold-mean coeffs for each measure at all
            % times (not just at start) 
            for j=2:thresh(i)
                t_k_mean_m(j,:)=t_k_mean_m(j-1,:)+exp(t_k_mean_raw(j,:));
            end  
            
            % remove the already assigned parameters from the constant list
            constants(1:thresh(i)*m_char_n(i))=[];
        
            % loadings for multi measure
            loading_m=loading(1:m_char_n(i));
            % remove the already assigned parameters from the loading list
            loading(1:m_char_n(i))=[];
        
            %%% simulated measure value for each individual and simulated factor
            %%% draw - excludes the constant term for multi-valued categorical
            %%% estimation
            %%% "f_ij" are simulated factor draws (common to all measures)
            %%% associated with a given factor
            sim_u=sim_factor*loading_m;


            % probability contribution of multi-valued measures: 
            % lower end
            sim_l_contrib_m=zeros(size(X_c_m,1),size(X_c_m,2));
            temp=bsxfun(@plus,-sim_u,t_k_mean_m(1,:));
            sim_l_contrib_m(X_c_m==0)=normcdf(temp(X_c_m==0));
            %thresholds in between
           for j=1:thresh(i)-1
                temp=bsxfun(@plus,-sim_u,t_k_mean_m(j,:));  
                temp1=bsxfun(@plus,-sim_u,t_k_mean_m(j+1,:));         
                sim_l_contrib_m(X_c_m==j)=normcdf(temp1(X_c_m==j))-normcdf(temp(X_c_m==j));
           end
           % upper end
           temp=bsxfun(@plus,-sim_u,t_k_mean_m(thresh(i),:));  
           sim_l_contrib_m(X_c_m==thresh(i))=1-normcdf(temp(X_c_m==thresh(i)));
           
           sim_l_contrib_m(isnan(X_c_m))=NaN;
            
           %%% Adding the multi-valued contributions to the overall factor likelihood 
           sim_l_contrib(:,b_char_n(i)+1:b_char_n(i)+m_char_n(i))=sim_l_contrib_m;
       
        end
        
        
        %%%%%%%%%%%% Simulated likelihood contributions of the CONTINUOUS
        %%%%%%%%%%%% measures (calculates sim_u and sim_l_contrib for those now)
        %nuisance parameters for the multi-valued measures associated with the
        %factor
        if cont_char_n(i) > 0
            
            %%% Gets the continuous measures for the factor 
            
            X_c_c=f_measures(:,b_char_n(i)+m_char_n(i)+1:b_char_n(i)+m_char_n(i)+cont_char_n(i));
            %nuisance parameters for the continuous measures associated with the
            %factor
            mean_c=constants(1:cont_char_n(i))';
            % remove the already assigned parameters from the constant list
            constants(1:cont_char_n(i))=[];

            % loadings for binary measure
            loading_c=loading(1:cont_char_n(i));
            % remove the already assigned parameters from the loading list
            loading(1:cont_char_n(i))=[];        

            % transforms the mean vector to a matrix so that it can be added to he
            % simulated loading contribution (the mean is common to all persons for
            % a given measure)
            mean_f=repmat(mean_c, rows,1) ;


            % simulated factor draws (common to all measures)
            sim_u=mean_f+sim_factor*loading_c;
            
            %%% "mu" just says that the error terms of each measure are assumed
            %%% mean zero (and sd sigma but that comes in later)
            mu=zeros(size(sim_u));
            %%% allows for a measure-specific error term 
            sigma=repmat(sd_cont',rows,1);
            sim_l_contrib(:,b_char_n(i)+m_char_n(i)+1:c_char_n(i)) = pdf('Normal',X_c_c-sim_u,mu,sigma);       
        
        end
        
    %%% cleans up variables for next iteration    
    f_measures(:,1:c_char_n(i))=[];
    sim_l_contrib_all_factors(1:rows,z:z-1+c_char_n(i))=sim_l_contrib;
    z=z+c_char_n(i);
    sim_factors(:,i)=sim_factor;
    clear alphas constants loading sd sd_cont X_c_b t_k_mean_b sim_u sim_l_contrib X_c_m t_k_mean_raw t_k_mean_m sim_l_contrib_m temp temp1 X_c_c mean_c mean_f mu sigma loading_b loading_m loading_c size_const
    
    end

   
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lottery Part %%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% All Types  %%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% parameter assignment for risk analysis
    
    theta1=theta(1);
    theta2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];        

    kappa1=theta(1);
    kappa2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];     

    sig_th1=theta(1);
    sig_th2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];  
    
 
    type1_3=NaN(3,n_type);
    p_type=zeros(n_type,1);
    
    sim_l_contrib_lottery=NaN(ind_n*draw_n,lottery_n);
    av_l_part=NaN(ind_n,n_type);
    theta_i=NaN(ind_n*draw_n,n_type);
    kappa_i=NaN(ind_n*draw_n,n_type);
    
    %%% weigh over unobserved types, each has its own set of intercepts for
    %%% the structural parameters of interest
    for tt=1:n_type
        %%% last type's proba is just 1-sum(other probas) so don't calculate
        if tt<n_type
            type1_3(:,tt)=theta(1:3);
            %%% probability needs to be between 0 and 1
            p_type(tt)=normcdf(theta(4));
            theta(1:4)=[];      
        elseif tt==n_type
            type1_3(:,tt)=theta(1:3);
            theta(1:3)=[];  
            p_type(tt)=1-sum(p_type);
        end
        
        %%% individual risk aversion based on sex and particular factor draws -
        %%% deterministic part - note that as above, the factor is determined
        %%% by the orthogonal draw AND by Xalpha - a function of observables

        %%% Even though the following parameters are individual-specific - thus
        %%% _i - they vary here within individuals due to the different factor
        %%% draws used in each row for the factor simulation
        theta_i(:,tt)=type1_3(1,tt)+sex*theta1+sim_factors*theta2_n;

        %%% individual trembling based on sex and particular factor draws -
        %%% normcdf() to bind it between 0 and 0.5 now (say that one can
        %%% make mistakes at most 50% of time, otherwise, probably
        %%% preferences misestimated and kappa takes on their role).
        kappa_i(:,tt)=0.5*normcdf(type1_3(2,tt)+sex*kappa1+sim_factors*kappa2_n);

        %%% individual variability of risk aversion based on sex and particular
        %%% factor draws 
        sig_th_i=sig_th_scale*normcdf(type1_3(3,tt)+sex*sig_th1+sim_factors*sig_th2_n)+sig_th_lower_limit;
        
        theta_i=min(maxima(1), theta_i);
        theta_i=max(minima(1), theta_i);
        kappa_i=min(maxima(2), kappa_i);
        kappa_i=max(minima(2), kappa_i);
        sig_th_i=min(maxima(3), sig_th_i);
        sig_th_i=max(minima(3), sig_th_i);        
        
        %%% compares an individual's level of risk aversion given each set
        %%% of simulated factor draws to the threshold level of risk
        %%% aversion proper to a given lottery choice task
        theta_xy_comp=kron(ones(ind_n*draw_n,1), theta_xy);
        theta_diff=(bsxfun(@plus,theta_xy_comp,-theta_i(:,tt)));
        temp=zeros(size(theta_diff));
        for i=1:lottery_n
            temp(:,i)=theta_diff(:,i)./sig_th_i;
        end
        %%% Probability of choosing the risky option i.e. probability that the
        %%% individual's risk aversion coeff is below the threshold for a
        %%% particular lottery. The difference in the coeffs is above
        %%% divided by the individual's sd of the error term so that the
        %%% normal cdf can be used here.
        %%% assumes normal errors on the risk aversion parameter:
        %%% theta~N(0,(sig_th)^2)
        p_risky=normcdf(temp);
        p_safe=1-normcdf(temp);

        % kappa percent of the time and individual makes the wrong choice.
        choice_risky=zeros(size(p_risky));
        for i=1:lottery_n
            choice_risky(:,i)=p_risky(:,i).*(1-kappa_i(:,tt))+p_safe(:,i).*kappa_i(:,tt);
        end
        %same as above and more efficient
        choice_safe=1-choice_risky;
        sim_l_contrib_lottery=(choice_risky.^(lottery_choice_ind)).*(choice_safe.^(1-lottery_choice_ind));
        
        
        
       %%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions from the various
       %%%%%%%%%%%%%%%%%%%%%% types of measures and lottery task choices

       sim_l_contrib_total=[sim_l_contrib_all_factors sim_l_contrib_lottery];


       % for numerical estimation purposes
       sim_l_contrib_total_plus=sim_l_contrib_total+1*10^(-20);
        
       
     
        % Individual likelihood contribution conditional on factor draws
        sim_ind_draw_contrib=prod(sim_l_contrib_total_plus,2);

        %%% Done so to efficiently store already estimated likelihood
        %%% contributions (in this case from factors and risk) by type.
        sim_l_contrib_1j(:,tt)=sim_ind_draw_contrib;
       
        % individual likelihood of measures and lottery choices averaged
        % over simulated factor draws
        av_l_contrib=cell2mat(accumarray(ind_id,sim_ind_draw_contrib,[],@(x){sum(x)}))/draw_n;
        
        % In order to avoid zeros for the log
        av_l_part(:,tt)=av_l_contrib+1*10^(-200); 
        
        sim_l_contrib_lottery=NaN(ind_n,lottery_n);
        clear sim_l_contrib_total sim_l_contrib_total_plus sim_ind_draw_contrib av_l_contrib sig_th_i theta_xy_comp theta_diff p_risky p_safe choice_risky choice_safe

    end

    %%%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions of the
    %%%%%%%%%%%%%%%%%%%%%%% various unobserved types and weigthing them by each type's
    %%%%%%%%%%%%%%%%%%%%%%% estimated population prevalence
    
    av_l=zeros(ind_n,1);
    for tt=1:n_type
        av_l=av_l+p_type(tt)*av_l_part(:,tt);
    end
    

    log_likelihood=-nanmean(log(av_l));
    
    %%% will capture the % of nan and inf in the final average likelihoods
    nan_pct=1-mean(isfinite(av_l));
    if nan_pct>.05
        log_likelihood=NaN;
    end
    
    %%% If at end the calculated proba of the last type is still <0,
    %%% discard estimate
    p_last=p_type(end);
    if p_last<0
        log_likelihood=NaN;
    end
    
end